REPORT  DOCUMENTATION  PAGE 


AFRL-SR-AR-TR-09-0138 


Public  reporting  burden  for  this  collection  of  information  is  estimeted  to  average  1  hour  per  response,  including  the  time  for  reviewing  inst  s 

dete  needed  end  completing  end  reviewing  this  collection  of  information  Send  comments  regarding  this  burdan  estimate  or  eny  othar  e;  j 

this  burden  to  Department  of  Defense,  Weshington  Headquerters  Services  Directorate  for  Informetion  Operations  and  Reports  (0704-01 


4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law  no  person  shell  be  subject  to  eny  penelty  for  felling  — , - «itly 

valid  QMS  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

03/20/09  Final 

3.  DATES  COVERED  (From  -  To) 

4/1/06  -  11/30/09 

4.  TITLE  AND  SUBTITLE 

Estimating  Neutral  Atmosphere  Drivers  using  a  Physical  Model 

5a.  CONTRACT  NUMBER 

FA9  550-06-1-0224 

5b.  GRANT  NUMBER 

FA9550- 06- 1- 022  4 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Tim  Fuller-Rowell,  Cliff  Minter,  Mihail  Codrescu,  Mariangel  Fedrizzi 

5d.  PROJECT  NUMBER 

1542571 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Colorado  at  Boulder 

572  UCB 

Boulder,  CO  80309-0572 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

1542571 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

DOD  AF  AFOSR 

875  North  Randolph  St. 

Arlington,  VA  22230-1768 

Dr.  Kent  Miller/NE 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


DISTRIBUTION  A:  Approved  for  Public  release 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

attached 


20090429225 


15.  SUBJECT  TERMS 

satellite  drag,  neutral  density,  physics  based  models,  data  assimilation 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Tim  Fuller-Rowell 

a.  REPORT 

public 

b.  ABSTRACT 

public 

c.  THIS  PAGE 

public 

8 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

303-497-5764 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z36.18 


Final  Report  March  2009 

Title:  “Estimating  the  Neutral  Atmosphere  Drivers  using  a  Physical 
Model” 

PI:  Tim  Fuller-Rowell 

Co-Investigators:  Cliff  Minter,  Mihail  Codrescu,  Mariangel  Fedrizzi 
Agency:  DOD  AF  AFOSR 
Award:  FA9550-06-1-0224 

In  the  initial  phase  of  this  study,  research  focused  on  implementing  a  physical  model  into 
an  ensemble  Kalman  filter  to  estimate  the  location  and  magnitude  of  the  upper 
atmospheric  heating.  Since  the  neutral  atmosphere  is  strongly  externally  driven  during 
geomagnetic  storms,  specifying  the  upper  atmospheric  heating  is  necessary  when 
describing  the  time-dependent  evolution  of  the  neutral  density.  Challenges  arise  when 
only  a  portion  of  the  upper  atmosphere  is  observable.  Because  the  upper  atmospheric 
dynamics  is  nonlinear  and  incompletely  observed,  determining  the  heating  location  from 
the  observed  atmospheric  response  is  not  straightforward,  as  shown  in  the  example 
Figure  1.  The  nonlinear  system  of  equations  can  provide  numerous  possible  solutions, 
only  one  of  which  is  correct.  Although  the  response  of  the  neutral  density  to  the  heating 
is  initially  linear  at  storm  onset,  the  system  becomes  increasingly  nonlinear  with  time  as 
illustrated  in  Figure  1.  The  true  heating  source  (blue  circles)  becomes  increasingly  more 
difficult  to  predict  (red  crosses)  as  the  storm  progresses. 


Figure  1 :  The  increasingly  nonlinear  response  to  the  upper  atmospheric  drivers  as  the 

geomagnetic  storm  progresses. 


The  physical  model  was  incorporated  into  the  Kalman  filter  system  in  two  phases  to 
better  tackle  the  challenges  in  dealing  with  the  nonlinear  heating/response  system.  The 
first  phase  used  a  simplified  model  that  numerically  solves  the  differential  equations  that 
describe  the  gravity  wave  propagation  of  the  upper  atmosphere  at  a  single  altitude.  Using 
the  simpler  model  to  solve  the  complex  nonlinear  dynamics  of  the  neutral  atmosphere 
provides  the  groundwork  for  building  the  full  system.  The  simpler  model  and  associated 
enKF  solution  allows  one  to  more  easily  tackle  the  challenges  in  estimating  a  nonlinear 
system.  The  capability  of  this  system  using  the  simple  model  is  illustrated  in  Figures  2 
below. 
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Figure  2:  A  simulation  of  an  unobserved  heating  source  (right)  and  the  enKF  solution  of 

the  unobserved  heating  (left). 


In  Figure  2,  a  heating  source  is  shown  by  a  simulated  density  spike  at  80  degrees  latitude 
and  1 80  degrees  longitude  in  the  left  panel  of  the  figure.  The  satellite  is  observing  a 
region  (the  outlined  white  satellite  track  covering  all  latitudes  at  360  degrees  in  the  left 
panel)  far  from  the  heating  source.  The  enKF  estimates  the  amount  of  heating  at  the 
unobserved  location  (at  1 80  degrees  longitude)  by  comparing  the  observed  response  (at 
360  degrees  longitude)  with  the  physical  model  prediction.  The  enKF  is  able  to 
iteratively  estimate  (right  panel)  what  the  heating  would  have  been  at  1 80  degrees 
longitude  based  on  the  observed  neutral  atmospheric  response  and  the  model  prediction. 
The  ability  to  solve  these  types  of  problems  is  important  because  they  resemble  the  actual 
complexity  of  the  system  as  well  as  the  inability  to  completely  observe  these  systems  due 
to  limited  satellite  coverage.  The  combination  of  the  nonlinear  dynamics  and  limited 
observation  coverage  makes  this  type  of  problem  very  difficult  to  solve.  The  system 
presented  here,  even  with  the  simplified  physical  model,  is  so  challenging  that  the  field  of 
nonlinear  optimization  and  other  related  fields  are  dedicated  to  solving  these  types  of 
problems.  The  ability  of  the  enKF  to  obtain  a  solution  for  this  nonlinear  system  indicates 


the  feasibility  of  the  solution  method  used  in  this  study  and  provides  a  framework  for  the 
full  physical  model  to  be  applied. 

The  second  phase  involved  replacing  the  simple  model  with  a  significantly  more 
advanced  numerical  model,  the  Coupled  Thermosphere-Ionosphere  Model,  CTIM 
(Fuller-Rowell,  et  al.  1996b).  CTIM  simulates  the  time-dependent  structure  of  the  wind 
vector,  temperature  and  density  of  the  neutral  thermosphere  by  numerically  solving  the 
non-linear  primitive  equations  of  momentum,  energy  and  continuity.  The  global 
atmosphere  is  divided  into  a  series  of  elements  in  geographic  latitude,  longitude  and 
pressure.  The  time  dependent  variables  of  southward  and  eastward  wind,  total  energy 
density,  and  neutral  temperature  are  evaluated  at  eaeh  grid  point  by  an  explicit  time¬ 
stepping  numerical  technique.  After  each  iteration,  the  vertical  wind  is  derived,  together 
with  temperature,  density,  and  heights  of  pressure  surfaces.  The  parameters  ean  be 
interpolated  to  fixed  heights  for  comparison  with  the  Challenging  Minisatellite  Payload 
(CHAMP)  or  the  Gravity  Recovery  and  Climate  Experiment  (GRACE)  density  data. 

The  implementation  of  the  physical  model  required  reformulation  to  operate  optimally  in 
a  recursive  data  assimilation  system  and  required  the  design  of  a  Kalman  smoother  to 
take  advantage  of  this  version  of  CTIM.  This  improved  system  overcame  the  major 
challenges  associated  with  solving  an  externally  foreed,  nonlinear  and  incompletely 
observed  neutral  density. 

Reformulating  CTIM  to  optimally  run  in  a  data  assimilation  system  has  been  a 
considerable  and  important  undertaking.  CTIM  was  originally  designed  for  investigative 
research  into  the  upper  atmospherie  dynamics  -  as  opposed  to  recursive  use  in  a  data 
assimilation  system.  Although  CTIM  is  regarded  as  an  advanced  model  in  the  field  of 
upper  atmospherie  research  and  has  proven  itself  for  operation  in  a  data  assimilation 
system,  its  forward-time-only  solution  inhibited  its  operation  in  a  recursive  system.  The 
data  assimilation’s  recursive  ability  is  a  necessity  for  solving  incompletely  observed, 
nonlinear  systems  where  the  external  forcing  is  not  entirely  known. 

This  recently  completed  recursive  data  assimilation  system  has  been  named  the  Global 
Recursive  Neutral  Density  Estimator  and  Loeator  (GReNDEL).  This  system  is  the  first 
of  its  kind  in  that  it  applies  an  advanced  physical  model  in  a  time-forward  and  -reverse 
recursive  system.  The  system’s  recursive  ability  accomplishes  several  goals  that  have  not 
been  possible  with  prior  systems,  whieh  includes:  (1)  a  complete  description  of  the 
location,  duration,  and  spatial  distribution  of  the  upper  atmospheric  heating  through 
observing  the  heating  response;  (2)  a  correct  nonlinear  description  of  the  upper 
atmosphere  response  onee  the  heating  has  been  specified  -  even  during  greatly  varying 
storm  conditions;  (3)  the  ability  to  calculate  conditions  in  unobserved  regions  through 
CTIM’s  nonlinear  dynamics;  and  (4)  the  ability  to  speeify  heating  anywhere  on  the  globe 
—  not  just  limited  to  the  eleetric  field  distribution.  This  system  is  a  major  step  in 
resolving  the  upper  atmospheric  heating  and  response  as  compared  with  traditional 
empirieal  model-based  systems,  whieh  must  rely  on  scalar  solar  and  geomagnetic  indices. 
The  results  of  this  progress  will  provide  an  improved  specification  of  the  neutral  density 
heating  distribution  beyond  the  previous  level  of  resolution. 


Solution  Method  of  GReNDEL 


Until  now,  data  assimilation  systems  have  been  limited  by  the  empirieal  model 
description  of  the  upper  atmospheric  nonlinear  dynamics  and  the  sealar  index  description 
of  the  varied  foreing.  Empirical  models  cannot  completely  describe  the  complex  chain  of 
events  that  connect  the  heating  and  the  complex  atmospheric  response  -  especially  if  the 
storm  continues  beyond  several  hours.  Furthermore,  the  use  of  sealar  geomagnetic 
indiees  to  describe  the  greatly  varying  heating  distributions  is  insufficient.  Small 
inconsistencies  in  the  heating  location  and  magnitude  ean  lead  to  vastly  different 
conclusions  since  the  upper  atmosphere  is  strongly,  externally  driven.  Furthermore,  these 
difficulties  are  exaeerbated  by  the  inability  to  observe  all  parts  of  the  globe  at  any  given 
time  with  many  of  the  new  single-satellite  in  situ  and  remote  observation  sourees.  Even 
global  ballistic  coefficient  estimation  does  not  have  a  sufficient  time  resolution  (several 
hours)  to  capture  rapidly  varying  storm  conditions  (minute-to-minute). 

The  recursive  physical  model  data  assimilation  system,  GReNDEL,  overcomes  the 
limitations  in  scalar  indices  and  empirical  modeling.  GRcNDL  improves  upon  the  solar 
and  geomagnetie  sealar  indices  by  instead  describing  the  location,  shape,  and  magnitude 
of  the  heating  source.  Once  this  more  complete  description  of  the  heating  is  obtained,  the 
physical  model,  CTIM,  numerically  solves  the  complex,  nonlinear  response  of  the  neutral 
atmosphere.  Recursive  iteration  of  this  system  provides  a  global  map  of  the  heating 
sources  -  even  in  regions  far  from  the  heating  souree  and  in  unobservable  regions. 

GReNDEL  accomplishes  these  tasks  by  recursively  solving  CTIM's  response  to  a 
particular  heating  source  (see  Figure  3).  Since  CTIM  can  now  operate  in  reverse-time, 
CTIM  ean  estimate  possible  heating  source  distributions  based  on  the  observed  response. 
Beeause  of  the  nonlinear  nature  of  the  system,  slight  differences  in  the  observed  response 
can  indicate  greatly  varying  possible  heating  conditions.  Because  of  the  sensitivity  and 
instability  of  the  nonlinear  system,  recursive  iterations  are  required  so  that  the  Kalman 
filter  ean  gradually  build  a  pieture  of  the  heating  location,  distribution  shape,  and 
magnitude.  The  result  is  a  pieture  that  is  much  more  descriptive  and  realistic  as 
eompared  with  scalar  indices  or  non-iterative  solution  methods. 

There  are  other  advantages  to  solving  the  system  in  this  recursive  manner.  Heating 
locations  at  other  locations,  beyond  the  observed  response  region,  can  be  solved.  This 
capability  is  neeessary  sinee  the  response,  due  to  the  propagation  of  gravity  waves,  is 
often  not  in  the  same  region  as  where  the  heating  occurred.  The  recursive  method  in 
GReNDEL  allows  one  to  accomplish  these  propagated  heating/responsc  solutions. 
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Figure  3:  Running  CTIM  in  reverse  time  and  recursively  correcting  heating/response 
to  estimate  density  as  well  as  heating  source  location,  distribution,  and  maginitude. 

CTIP  Validation  against  Multiple  Data  Sources 

Having  established  a  data  assimilation  infrastructure  using  CTIM,  the  next  phase  was  to 
improve  the  physical  model  by  incorporating  a  more  realistic  global  ionosphere- 
plasmasphere,  remove  biases,  and  improve  the  code  to  follow  the  response  and  recovery 
time-scales  to  geomagnetic  activity.  CTIM  used  a  simplified  ionosphere  limiting  the  use 
of  plasma  observables  within  the  data  assimilation  scheme.  Detailed  comparison  of  the 
more  advanced  code  that  includes  a  global  ionosphere  plasmaspherc,  the  Coupled 
Thermosphere  Ionosphere  Plasmasphere  electrodynamics  model  (CTIPe),  enable  the 
physical  model  to  be  validated  against  neutral  and  plasma  data  sources.  In  the  current 
tests,  the  data  sources  include  the  CHAMP  neutral  density  data,  the  European  incoherent 
scatter  radar  facility  (EISCAT)  ionospheric  data,  and  total  electron  content  from  the 
operational  US-TEC  data  assimilation  model.  Detailed  comparison  enabled  design  of 
some  significant  improvements  in  the  code  to  remove  model/data  biases  and  model  time 
constants. 

Figure  4  shows  a  comparison  of  CTIPe  with  eight  days  of  a  year-long  operation  of  the 
EISCAT  facility.  The  model  is  able  to  follow  the  day-to-day  changes  during  April  2005 
reasonably  well  with  minimal  biases.  However,  the  model  predicted  a  winter  anomaly  at 
high  latitudes,  which  was  not  observed  by  EISCAT,  which  introduced  a  seasonal  bias. 
This  bias  was  removed  by  introducing  more  realistic  seasonal/latitude  dependence  in 
solar  dissociation  rates. 

Detailed  comparison  of  CTIPe  with  CHAMP  during  six  months  in  2005  also  revealed 
that  the  recovery  to  geomagnetic  storms  was  too  slow  in  the  numerical  simulations.  The 
physical  process  affecting  the  recovery  times  is  largely  due  to  the  radiative  cooling  by 
nitric  oxide.  If  the  recovery  rate  is  too  slow  the  magnetospheric  energy  input  required  to 
maintain  the  observed  density  was  also  correspondingly  too  small.  The  NO  cooling  rate 
during  a  geomagnetic  storm  is  controlled  by  too  processes.  Firstly,  the  increase  in 
temperature  in  the  presence  of  the  background  solar  produced  NO  increases  the  cooling 
rate.  Secondly,  the  NO  production  at  mid  to  high  latitudes  during  a  geomagnetic  storm  is 


controlled  by  auroral  electron  precipitation,  which  further  increases  the  cooling  rate.  The 
latter  process  was  underestimated  in  CTIPe  and  in  CTIM. 
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Figure  4.  Comparison  of  EISCAT  plasma  density  (black)  with  a  physics-based  model 

To  capture  the  recovery  times  more  accurately  an  empirical  NO  model  from  Marsh  et  al., 
based  on  SNOE  observations,  was  included  as  a  module  in  order  to  simulate  the  observed 
height/latitude  storm  time  changes  in  NO  more  realistically.  The  model  is  able  to  improve 
the  storm-time  changes  in  NO  and  radiative  cooling  during  the  response  and  recovery 
phase  of  a  geomagnetic  storm.  The  added  advantage  is  that  the  empirical  model  is 
computationally  fast  and  robust  under  all  forcing  conditions. 

Figure  5  shows  a  comparison  of  CTIPe  neutral  mass  density  with  the  CHAMP  satellite 
data  for  the  first  fifteen  days  of  2005.  Model  biases  are  not  present  for  this  period  and  the 
response  and  recovery  time  scales  are  realistic.  This  more  accurate  physics-based  model 
can  be  used  as  the  forward  model  in  the  data-assimilation  scheme  and  for  subsequently 
forecasting  the  system.  The  specific  state  of  the  thermosphere-ionosphere  is  dependent  on 
the  energy  sources:  particularly  solar  radiation  and  magnctosphcric  Joule  and  particle 
heating,  together  with  internal  mixing  processes  from  turbulence  and  global  circulation 
(the  thermospheric  “spoon”). 

Extensive  comparison  of  CTIPe  has  also  been  done  with  the  total  electron  content  (TEC) 
over  the  CONUS  from  the  US-TEC  operational  model.  Validation  of  US-TEC  itself  has 
revealed  an  accuracy  of  less  than  2  TEC  units  for  vertical  TEC,  so  providing  an  excellent 
data  source  to  validate  the  physics  based  model.  CTIPe  has  been  running  in  a  real-time 
mode  for  about  a  year  to  enable  comprehensive  comparison  with  the  data  to  establish  the 
baseline  accuracy  through  the  seasons. 

The  structure  of  the  code  was  also  modularized  and  the  code  was  made  more  robust  to 
enable  stronger  forcing  terms  to  be  applied.  With  the  increase  in  radiative  cooling  it  was 
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necessary  to  include  much  larger  geomagnetic  sources,  particularly  Joule  heating,  in 
order  to  follow  the  time  history  of  neutral  density  during  a  geomagnetic  event. 
Modularizing  the  code  enable  the  physical  processes  in  the  more  sophisticated  CTIPc 
code  to  be  ported  to  CTIM  if  required,  which  is  a  simpler  faster  code  suitable  for  data 
assimilation.  In  addition,  modules  that  are  not  required  to  be  self-consistent  can  be 
replaced  by  empirical  algorithms.  Validating  and  estimating  drives  for  one  parameter  in  a 
physical  model  (e.g.  Nc)  provides  a  means  to  specify  and  forecast  another  parameter, 
such  as  neutral  density  for  satellite  drag. 
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Figure  5.  Comparison  of  CHAMP  satellite  orbit-averaged  neutral  density  (black)  with  a 

physics-based  prediction. 


The  next  step,  which  will  be  continued  under  the  new  NADIR  MURI  initiative,  will  be  to 
incorporate  the  more  sophisticated  and  validated  physical  model  into  the  ensemble 
Kalman  filter  data-assimilation  scheme. 
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